In this tutorial we will learn how to
The MITK-MSI software provides a wrapper to the popular MCML approach to simulate how light travels through tissue. This wrapper can be found in mc/sim.py. In this tutorial we will utilize our tissue model which uses this wrapper to create the reflectance spectra.
As a prerequisit, you need a MCML monte carlo simulation which uses the format specified here. I tested this software with the GPU accelerated version which can be found here.
In [1]:
# 1.1 create spectra - setup simulation environment
# some necessary imports
import logging
import numpy as np
import os
# everything related to the simulation wrapper
from mc import sim
# the factories create batches (tissue samples) and suited tissue models
from mc import factories
# function which runs simulations for each wavelength
from mc.create_spectrum import create_spectrum
# Where does your monte carlo simulation executable resides in?
MCML_EXECUTABLE = "/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/gpumcml.sm_20"
# The MCML needs a simulation input file, where shall it be created?
MCI_FILENAME = "./temp.mci"
# filename of the file with the simulation results. Due to a bug in GPUMCML will always
# be created in the same folder as the MCML executable
MCO_FILENAME = "temp.mco"
# The wavelengths for which we want to run our simulation
WAVELENGTHS = np.arange(450, 720, 2) * 10 ** -9
# we want to create standard colonic tissue as specified in the IPCAI 2016 publication
# "Robust Near Real-Time Estimation of Physiological Parameters from Megapixel
# Multispectral Images with Inverse Monte Carlo and Random Forest Regression"
factory = factories.ColonMuscleMeanScatteringFactory()
# if you want to create data from the generic tissue mentioned in the paper, choose:
#factory = factories.GenericMeanScatteringFactory()
# create a simulation wrapper
# the simulation wrapper wraps the mcml executable in python code
sim_wrapper = sim.SimWrapper()
# our simulation needs to know where the input file for the simulation
# shall resign (will be automatically created)
sim_wrapper.set_mci_filename(MCI_FILENAME)
# also it needs to know where the simulation executable shall lie in
sim_wrapper.set_mcml_executable(MCML_EXECUTABLE)
# create the tissue model
# it is responsible for writing the simulation input file
tissue_model = factory.create_tissue_model()
# tell it where the input file shall lie in
tissue_model.set_mci_filename(sim_wrapper.mci_filename)
# also set the output filename
tissue_model.set_mco_filename(MCO_FILENAME)
# tell it how much photons shall be simulated. Will be set to 10**6 by standard,
# this is just an example
tissue_model.set_nr_photons(10**6)
In [2]:
# 1.2 create spectra - create tissue samples for simulation
# setup batch with tissue instances which should be simulated
batch = factory.create_batch_to_simulate()
# we want to simulate ten tissue instances in this example
nr_samples = 10
df = batch.create_parameters(10)
# lets have a look at the dataframe. Each row corresponds to one tissue instance,
# each tissue instance is defined by various layers, which all have certain parameters
# like e.g. oxygenation (here sao2)
df
Out[2]:
In [3]:
# 1.3 create spectra - run simulation
# add reflectance column to dataframe
for w in WAVELENGTHS:
df["reflectances", w] = np.NAN # the reflectances have not been calculated yet, thus set no nan
# for each instance in our batch
for i in range(df.shape[0]):
# set the desired element in the dataframe to be simulated
tissue_model.set_dataframe_row(df.loc[i, :])
logging.info("running simulation " + str(i) + " for\n" +
str(tissue_model))
reflectances = create_spectrum(tissue_model, sim_wrapper, WAVELENGTHS)
# store in dataframe
for r, w in zip(reflectances, WAVELENGTHS):
df["reflectances", w][i] = r
# clean up temporarily created files
os.remove(MCI_FILENAME)
created_mco_file = os.path.join(os.path.split(MCML_EXECUTABLE)[0], MCO_FILENAME)
if os.path.isfile(created_mco_file):
os.remove(created_mco_file)
# Hooray, finished,
# now our dataframe also contains reflectances for each tissue instance:
df["reflectances"]
Out[3]:
In [4]:
# 2.1 analyse spectra - plot reflectances
# the usual settings for plotting in ipython notebooks
import matplotlib.pylab as plt
%matplotlib inline
# let's have a look at our reflectances
df["reflectances"].T.plot(kind="line")
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
In [5]:
# 2.1 analyse spectra - show distribution of blood volume fraction (vhb) and sao2
# now we need some special pandas functions
import pandas as pd
# we're interested in the distribution of vhb and sao2 in the first layer (layer0)
df_vhb_sao2 = df["layer0"][["vhb", "sao2"]]
# plot a scatter matrix showing the distribution of vhb and sao2.
# of course, with this little data this does not really make sense,
# however it is a useful tool for analysis if much data is available
pd.tools.plotting.scatter_matrix(df_vhb_sao2, alpha=0.75, figsize=(6, 6))
plt.show()
In [6]:
# 3.1 manipulate spectra - apply sliding average
# in 3.1 and 3.2 we will adapt the generated spectra to an imaginary imaging system
# This system has filters with 20nm bandwith (taken care of in 3.1)
# and takes multispectral images in 10nm steps (taken care of in 3.2)
# the module mc.dfmanipulations was written to provide some basic,
# often needed manipulations of the calculated spectra
# all dmfmanipulations are performed inplace, however, the df is also returned.
import mc.dfmanipulations as dfmani
# first copy to not lose our original data
df2 = df.copy()
# We apply a sliding average to our data. This is usefull if
# we want to see e.g. how the reflectance was recorded by bands with a certain width
# a sliding average of 11 will take the five left and five right of the current reflectance
# and average. Because we take 2nm steps of reflectance in our simulation this means
# a 20nm window.
dfmani.fold_by_sliding_average(df2, 11)
# lets again plot the reflectances
df2["reflectances"].T.plot(kind="line")
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
# we can see that the bump at 560nm is "smoother"
In [7]:
# 3.2 manipulate spectra - select certain wavelenghts
# our imaginary imaging system takes images in 10nm steps from 470 to 660nm
imaging_system_wavelengths = np.arange(470, 670, 10) * 10**-9
df3 = df2.copy()
dfmani.interpolate_wavelengths(df3, imaging_system_wavelengths)
# let's look at the newly created reflectances
df3["reflectances"].T.plot(kind="line", marker='o')
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
In [12]:
# that's it, folks. If you want, you can save the created dataframe easily to csv:
df.to_csv("results.csv", index=False)
In [ ]:
In [ ]: